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We show that for two identical neuronal oscillators with strictly postive phase restting curve, 
isochronous synchrony is an unstable attractor and arbitrarily weak noise can destroy entraiment 
and generate intemittent phase slips. Small inhomogeneity-mismatch in the intrinsic firing rate of 
the neurons- can stabilize the phase locking and lead to more precise relative spike timing of the two 
neurons. The results can explain how for a class of neuronal models, including leaky itegrate-fire 
model, inhomogeneity can increase correlation of spike trains when the neurons are synaptically 
connected. 
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Synchronization observed frequently in the vast variety 
of physical, chemical, industrial, and biological complex 
systems, from coupled pendulum clocks to neuronal pop¬ 
ulations in the nervous system [1-10]. In these systems, 
ability to exhibit synchronous or phase locked oscilla¬ 
tions, is the foundation of the emergent behaviors which 
is the basis for the functionality of the system. While the 
existence of robust synchronization is important in real 
systems with different sources of noise and uncertainties 
in parameters, stability of this behavior has a central 
importance [4]. Recordings of multi-neuron spike trains 
have revealed significant interdependencies between the 
firing of different neurons in a population [11-16]. Syn¬ 
chronous oscillations are found in many brain regions and 
excessive synchrony is a hallmark of neurological disor¬ 
ders such as epilepsy and Parkinson’s disease [17]. Func¬ 
tional role of the correlation in neural coding has been 
debated in recent years [18-23]. Synchrony itself may 
encode information directly [10, 14, 16, 24-27]. Syn¬ 
chronous firing of the neurons in one region serves to re¬ 
liably transmit signals to upstream regions [28-30] , while 
synchrony between different regions can prepare dynamic 
channels for communication [31-33] and also undelies fea¬ 
ture binding [34]. Beyond the functional role, it is also 
important to understand how correlation and synchrony 
depend on biophysical parameters of the neurons and the 
network. 

Correlation between spike trains of neurons can arise 
from shared input they receive from other neurons [35- 
39], or from presence of direct synaptic connections be¬ 
tween neurons [40-43]. In both cases the collective state 
of the system depends on the parameter of neurons, 
e.g., firing rate and the type of the excitability of neu¬ 
rons [44], and the parameters of the connections such 
as delay [45]. Physiological heterogeneity can destabilize 
both coupling-induced and correlation-induced synchro¬ 
nization [46-48]. In the classical models of synchroniza¬ 
tion, collective state of a system of coupled oscillators is 
determined by outcome of rivalry between synchronizing 
effect of connections and desynchronizing effect of inho¬ 


mogeneity [3], but there are examples of the systems in 
which synchrony is enhanced by inhomogeneity [49, 50]. 
Recently we have shown that small inhomogeneity can 
increase correlation between spike trains of two coupled 
neurons [51]. In this sudy we give a general framework for 
the correlation of coupled phase oscillators with a given 
phase sensitivity. We show that for identical pulse cou¬ 
pled type-I oscillators, synchronized state is an unstable 
attractor and arbitrarily weak noise can destabilize this 
state and the spiking of two neurons exhibit intermittent 
phase slips between epochs of locking. Small inhomo¬ 
geneity in firing rates can stabilize the system by provid¬ 
ing an asymmetric basin of attraction around the stable 
phase-locked state. This in turn results in a sharper PDF 
for the time difference between spikes of the two neu¬ 
rons in presence of noise. We have also shown that while 
for the model neurons with biologically realistic phase 
response curve (PRC), the time difference between the 
spikes of two neurons in the stable state increases with 
inhomogeneity, in the case of LIF neurons, they lock in 
almost zero phase lag for sufficiently small values of in¬ 
homogeneity. By solving Fokker-Planck equation we also 
find the most probable phase difference between spike 
times of the two neurons and will show that it does not 
coincide with the stable point of the deterministic equa¬ 
tions. 

Our model comprises two bidirectionally coupld neu¬ 
ronal oscillators recieving suprathreshold constant cur¬ 
rents (/i and I 2 with mismatch A/) as well as indepen¬ 
dent stochastic inputs. The evolution of the state vector 
of the oscillators Xi,i = 1,2 can be descibed by 

Xi = F{Xi) cgi2Gi2{Xi,X2) -I- /i -I- cr^i(t) (1) 

X 2 = F{X2) + eg2iG2i{X2, Xi) + I 2 + o'^2{t), (2) 

where F governs the internal dynamics of the neurons, G 
detemines the synaptic connections, ^ is Gaussian white 
noise with zero mean and unit variance, and e and cr are 
small values which scale strength of the couplings and the 
stochatic inputs, respectively. We assume the each of the 
unperturbed systems Xi = F{Xi) has an asymptotically 
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stable limit cycle, Xo(t) = X^it + T), so that a phase 
variable can be defined in vicitnity of the limit cycle. In 
the regime of weak coupling and weak noise we can apply 
the standard phase reduction [52-54] to the Langevin 
equations above. The system is then can be described by 
a set of Ito stochastic differential equations: 

9i = LOi + egi 2 Z { 6 i) G { 61 , 92 ) 

+ /DiZ(0i)Ci(t) 

62=102 + sg 2 iZ { 62 ) G { 62 , 61 ) 

+ V^Z{ 62 )^ 2 {t) (3) 


where Z{ 6 ) is the infinitesimal phase-response curve 
(PRC) [55]. We assume that the natural frequencies 
have a small difference coi — L 02 = Aw and the noise 
and coupling influence only the first (voltage) variable 
of the state vector of the neural oscillators. In our 
model the neurons communicate via pulsatile signals 
Gij = '^^S{t — tj-), where <5 is Dirac’s delta function 
and is the instant of firing of the neuron j. These 
pulses idealize the communcation signals which are short 
compared to the intrinsic time scale of the oscillators and 
are used to model diverse systems such as populations 
of flashing fireflies and plate tectonics in earthquakes, 
as well as networks of spiking neurons in the brain [56- 
61]. It is assumed that the mismatch, coupling and noise 
terms are of the same order, sufficiently weak such that 
the intrinsic dynamics of isolated identical phase oscilla¬ 
tors is dominant. 

Using the method of averaging [62] we derive the the 
equation of motion for the phase difference (j) = 61 — 62 ' 

^ = e[Auj - gi 2 H {(j}) + g 2 iH {-(f))] 

+ a^y/2De g {t) (4) 


where H {(j)) = Z Z (J) G {t,i + </>) di and the coeffi- 

cient = ( ^ /o [^ (^] comes from averaging 

the noisy phase equations [52] . Here g {t) = 

is itself a Gausian white noise with zero mean and unit 

variance. 

We restrict the study to type-I oscillators and first dis¬ 
cuss on the deterministic version of Eq. 4 with D = 0. 
If H is an even function of (p, e.g. for QIF oscilla¬ 
tors, the effective coupling term would be AgH{<j)) with 
Ag = 521 — gi 2 - In this case the most effective connection 
is a unidirectional one and the symmetric connection has 
no effect on the relative dynamics of the oscillators. Note 
that for the oscillators with an oblique PRC, e.g. the LIF 
oscillators, the coupling term can be non-zero for sym¬ 
metric connections (see suplementary material Fig. SI). 

The fixed point of Eq. 4 with D = 0 is the solu¬ 
tion of Aw = gi 2 H {(p) — g 2 iH {—(p). For QIF oscilla¬ 
tors Z {p) = I — cos {p) with asymmetric connections 
Ag 7 ^ 0, when the oscillators are identical Aw = 0, the 






FIG. 1. (A,B) Representative examples of the evolution of 
the phase difference of two neurons for three different values 
of mismatch in intrinsic frequencies. Larger values of mis¬ 
match have led to fewer phase slips. In A neurons are phase 
oscillators with canonical type-I phase sensitivity and in B the 
results are presented for LIF neurons. (C) The mean scape 
time is plotted against frequency mismatch. Increasing ef¬ 
fective coupling constant Ag = gi — g 2 the maximum scape 
time is seen in largar values of frequency mismatch. In (D) 
the ratio of the firing rates of the coupled neurons is plotted. 
For large values of mismatch the fixed point of 1 : 1 locking 
vanishes. 


zero-lag synchrony </> = 0 is an unstable attractor and 
in the absence of noise, the oscillators can synchronize 
isochronously. But a waek noise can destroy synchrony 
and lead to phase slips. Mismatch in the intrinsic firing 
rates of the neurons, stabilizes the fixed point through a 
saddle-node bifurcation while moves the fixed point away 
from zero. For small mismatch, this provides an asym- 
matric basin of attraction which is vulnerable to suffi¬ 
ciently large perturbations in one direction around the 
fixed point. In the presence of noise the system shows 
epochs of intermittent locking between which the rela¬ 
tive phase of the oscillators slips by one cycle, while the 
mean scape time from locked states increases with fre¬ 
quency mismatch (see Figs. lA and B). The maximum 
mean scape time from the locked state, is seen in a certain 
value of mismatch (Fig. IC) and for larger mismatches 
Aw > Max{gx 2 H {p) — g 2 iH (—<(>)}, the fixed point cor¬ 
responding to 1 : 1 locked state will disappear through 
another saddle-node bifuracation (Fig. 3). To give more 
concrete results on the impact of the inhomogeneity on 
the correlation of the spike trains of the neuronal oscil¬ 
lators in presence of noise, we derive the Fokker-Planck 
equation for the distribution of the phase difference of 
two neurons, described by Eq. 4. We rewrite Eq. 4 in a 
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more closed form 

^ {(j)) + a^V2D£T] (t) (5) 


where T (</>) = ^ 


T 


4 COS(j) 


The corresponding 


Fokker-Planck equation takes the form: 


( 6 ) 


where p (</), t) is the distribution of the phase differences. 
Stationary solution of this equation with periodic bound¬ 
ary condition is: 



where M {(j)) = ^ /q^ T ( ^) d^- Also, N is the normaliza¬ 


tion factor so that p {(j)) d(j) = 1 and a = is the 
ratio of noise intensity to the coupling strength [63]. 

Figure 2A shows the steady state phase difference dis¬ 
tribution for different values of the frequency mismatch 
for QIF neuronal oscillators. It can be seen that the 
distribution becomes narrower (with a more pronounced 
peak) with increasing frequency mismatch while the neu¬ 
rons remain in 1 : 1 locked state, i.e. for the mismatch in 
the range 0 < ^ ^. This reflects a larger basin of at¬ 

traction for the focked state when mismatch is increased 
from zero. Furthermore, the asymmetry of the basin of 
attraction causes the distribution of the phase differences 
not to peak in the fixed point of the deterministic equa¬ 
tion, determined by F (^*) = 0. In turn, in presence of 
noise the location of maximum phase difference satisfies. 


F (<(>*)= a 




( 8 ) 


which is derived by taking the derivative of p (</>) with 
respect to (p, equal to zero. The location of the most 
probable phase difference as a function of mismatch, de- 
trermined by Eq. 8, is plotted in Fig. 2C for different 
values of noise amplitude as well as for noiseless system 
which shows the location of the fixed point. Presence of 
noise inclines the distribution to larger phase differences 
for small values of frequency mismatch. The maximum 
difference between the location of most probable phase 
difference between noiseless state and noisy state is seen 
near ^ = 0 or which reflects the most asymmetric 
basin of attraction for the locked state and in turn the 
locations coincide when p* = 7r/2 where the basin of at¬ 
traction is symmetric. 



FIG. 2. (A) The steady state phase difference distributions 
p{4>) for three levels of heterogeneity. Distributions have be¬ 
came narrower as mismatch is increased. Solid lines show 
the analytic result Eq. S28 and the bar graph presents the 
numerical results by direct integration of Eqs. S7. Dashed 
vertical lines show the position of the fixed points of deter¬ 
ministic equations. (B) The maximum value of p (p) is plotted 
against frequency mismatch for two different values of the the 

D(t^ 

ratio of noise strength to effective coupling a = (C) Tke 

most probable phase difference are shown for two values of a. 
They don’t concide with the location of fixed point of the 
deterministic equations (black curve). 


Pfeuty et al. (2005) have introduced a variable Si (t) 
which is equal to I /S when a neuron has fired a spike in a 
time bin of size S about time t and is equal to 0 otherwise 
[64]. For sufficiently small S the time average of Si is 
the average firing rate of neuron i. It is shown that the 
normalized cross-corellogram (CC) of this variable which 
is the density of probability that neuron 2 to fire a spike 
in a time bin of size S a delay r after a spike of neuron 1, 
is related to the phase difference probability distribution 
function p {(p) through 



{Si {t) S 2 {t -I- t)) 
{Si{t)){S2{.t)) 


(9) 


where (...) indicates averaging over time. A peak in CC 
at a time lag r shows phase locking of the activity of the 
neurons. The sharper CC is indicator of a tighter locking. 
To illustrate this effect we provide an expression for the 
maximum value of the distribution fuction (or CC) as a 
function of frequency mismatch: 


Cm 


ax 



= a 


r 


1 — e 





( 10 ) 
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FIG. 3. The cross-correlogram of spike trains of two LIF 
neurons C (r) shows that in presence of the mismatch cross 
correlation is increased. The level of maximum correlation is 
shown in the inset to highlight the increase due to the inho¬ 
mogeneity. 


The above equation is the same as Eq. 8 by substituting 
p {(jf) with C • Figure 2B shows maximum value of 
cross-correlation versus frequency mismatch for different 
values of noise to coupling ratio which is resulted from 
direct integration of Eq. 8. The result shows that the 
maximum cross-correlation of the spike trains of the os¬ 
cillators would be also maximum when the neurons are 
not identical. This is a consequence of more precise rela¬ 
tive spike timing of the two neurons in presence of inho¬ 
mogeneity. 

In this study we have shown that for two synapti- 
cally connected neuronal oscillators, more precise relative 
spike timing can be achieved when the neurons receive 
different levels of inputs and have different intrinsic fir¬ 
ing rates. Consequently, cross-correlation of spike trains 
of the neurons increases in presence of mismatch in in¬ 
trinsic firing rates of neurons. While the results are pre¬ 
sented for neuronal oscillators, they can find application 
in general context of coupled limit cycle oscillators. 
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Supplementary Material for 
“Stabilizing synchrony with heterogeneity” 


WEAKLY COUPLED OSCILLATORS 

Our model comprises two bidirectional coupled neurons receiving suprathreshold constant currents as well as un¬ 
correlated stochastic inputs. The general form of equations describing this model is given by 

Xi = F{Xi) J- egi2Gi2{Xi, X 2 ) -I- /i + cr^i (t) 

X 2 = F{X2) + eg2iG2i{X2, Xi) 12 + cr^2{t) (SI) 

where Xi is a N-dimensional state vector containing the membrane potential and gating variables. For example in 
the multicompartmental Hodgkin-Huxley (HH) model, X = [V, m, h, n]^ and in the single compartmental Leaky 
Integrate-and-Fire (LIF) model, X = V. F {X) defines the internal dynamics of the neuron i. Gtj {Xi,Xj) determins 
functional form of synaptic connection from neuron j to neuron i. For example in the case of pulse coupled synapses, 
as we used in this letter it would be 

G.,[x,{t),X,{t)) (S2) 

n 

where t" is the instant of spike of neuron j and S (t) is the Dirac’s delta function, gij shows the synaptic weight 
from neuron j to neuron i which is scaled by a small factor e so that the weakly coupled oscillators approximation is 
valid. 

Each neuron receives suprathreshold constant current li with mismatch A/ = Ii — I 2 as well as an independent 
Gaussian white noise characterized by its mean, auto- and cross-correlation with the other neuron’s input 

(C. {t)) = 0, (S3) 

= (S4) 

We assume that at the absence of synaptic connections and noise, dynamics of isolated neurons given by 
A, = F{X,) + h 


(S5) 
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have a T-periodic limit cycle solution Xq (t). We assumed that the magnitude of the coupling term scaled with a 
small coefficient e and also the amplitude of the noise is such that the variance of the noise and the strength of the 
coupling is in the same order, so tr = V De, where D = O {1). In such a weak coupling and weak noise regime the 
dynamics of the neurons can be approximated by defining a single phase variable around limit cycle. So we define a 
phase variable, 6 i (t) in the vicinity of unperturbed limit cycle for each oscillator and reduce high dimensional Eqs. 
SI to two scalar equations for the evolution of the phase 


01 = CCi + £512^ (0i) G (01,02) + (0i) a (t) 

02 = W2 + Sg2lZ (02) G (02, 0l) + VDeZ (02) ^2 {t) 


(56) 

(57) 


where Z (0) is the infinitesimal phase response curve (PRC). We assumed that neurons natural frequencies have small 
difference, eAu! = W 2 — wi, where e confirms that the mismatch is of the order O (e). Assuming Xi (t) ~ Xq {Oi (t)) 
and by the change of variable 6 i (t) = uit + 4'i (t), the equations for the evolution of the relative phase of the oscillators 
read: 


(j)i = - 
02 = + 


eAoj 

2 

eAuj 


+ egi2Z (0i) G (01,02) + y/lXZ (0i) ^i (t) 
+ eg2iZ (02) G (02,0i) + ViXZ (02) ^ (t). 


We exploit the fact that e is small to further reduce Eqs. S9. With a system of the form 
i = ef {x,t). 

Averaging theory states that in Eq. SIO, x (t) can be replaced by its average over a period x and 


1 

x = e— J f{x,t)dt. 


By applying averaging method on the Eqs. S9 we have 
Aw 

T 

Aw 


01 = + £^2" (02 - 0i) + cr^v^^i (f) 


02 = +£- 


+ £ 


(01 - 02) + a^'/D£^2 (t) 


(58) 

(59) 

(510) 

(511) 

(512) 

(513) 


where the term ^ /q [-^ (^] originates from averaging the noisy phase equations, and Z {(j)i — (j)j) 

comes from 


H (02 - 0l) = 




^(0i)G(0i,02)d0i 


Z(0i)5(02)d0i 


Z {t + 4ii) S {t + (j>2) dt 

^ Jo 

I F ~ 

- J Z (t) (5 (< - 01 + 02) di 


Z (02 - 0l) 
T 


(514) 

(515) 

(516) 

(517) 

(518) 


without loss of generality we assumed the phase is normalized so that O<0<T,i.e.,w = l. By defining 0 = 02 — 0i, 
we derive the following equation for the phase difference 


0 = e(^Auj + g 2 iH {-(/)) - gi 2 H{(j))j + a 2D eg (t) 


(S19) 


where {t) = ^2 — and g{t) is a Gaussian white noise with zero mean and unit variance. 

The have used two model neurons in ot study: Canoncal type-I oscillators with y 2’(0) = 1 — cos(0) and LIE 
oscillators which is described by u = / — u for t < VTh{= 1) and limr^o+v{t + r) = 0 with 


Z{(j)) = yexp( — \ 0 < 0 < 27r, 

1 \uj / 


(S20) 






3 


F(0) = As (1 - COS0) F(0) = ^ [(s - As) - (s + As)e 



FIG. SI. Examples oiT x F function, F {(/}) = gi 2 H {4>) — g 2 iH for QIF oscillator (left) and LIF oscillator (right) for 

different values of difference of the synaptic strengths. For the LIF oscillators with an uneven PRC the coupling term can be 
non-zero for symmetric connections whereas for QIF oscillators with an even PRC the effective coupling is determined by the 
difference of two reciprocal synaptic constants. 


where w = 2Tr/[logi — log{I — I)]. First we focus on the deterministic case of Eq. S19 with D = 0. For QIF oscillator, 
Z {(j)) = I — cos {(j)) is an even function of (j). Therefore it reduces to 


= e 


Aw' 


T ' 


1 — cos 


(S21) 


where Ag = g 2 i—gi 2 is the effective coupling constant. In this case the most effective coupling is that which maximizes 
Ag. i.e., a unidirectional one and the symmetric connection leads to zero effective coupling Ag = 0. But for LIF 
neurons with uneven PRC Eq. S19 takes the form 


= e 


A , ^921 

Aw-b -j^exp 
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0 < d) < 27r. 


(S22) 


Note that for the oscillators with an oblique PRC, e.g. the LIF oscillators, the effective coupling term can be non-zero 
for symmetric connections (see Fig. SI). The fix points of general equation S19 with D = 0 are the cross points of 
horizontal line TAuj and the curve described by F{(j)) = gi 2 Z{(l)) — g 2 iZ{—(j)). 

For type-I phase oscillators we rewrite Langevin Eq. S19 as 


^ = eZgV {(j)) -b a^\/2Der] (t) 


(S23) 
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with F (cj)) = 
is 

^ {cj), t) = -e^ [F {^) p {^, t)] + (cjlDe) (0, t). 

The stationary phase difference distribution satisfies 
dpo (</>, t) 


. Corresponding Fokker-Planck equation for the phase difference distribution p{cj), t) 


dt 


= 0 , 


with the solution 
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where 
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(S24) 
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(S27) 
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iV is a normalization factor so that f„ p{cj))d(l) = 1, and a = 


jg fy ~ A fS-tio of noisc intensity to the coupling 

strength. The constant A can be determined by the periodicity condition of po, that is, po (0) = po (T). Therefore 
the final form of stationary solution is 
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- Jo e 


T -M 




Jo 


+ 1 


(S28) 


In Figure 2A we have shown the result of analytic solution for steady-stat phase difference distributions S28 and 
that of direct numerical integration of of phase differential equations S7. For solving Eq. S28, we have used double int 
function of MATLAB. In simulation, we integrate Eqs. S7 with Euler method and save spike times of each neuron. 
Then we have used hist function in MATLAB to plot p. 

It has been shown that in the weak coupling and weak noise limit, the cross-crologram (CC) and the phase difference 
probability distribution function, p((^), are related by 


CM = p(J). 


(S29) 


The most probable phase difference of spiking of the neurons (location of the peak of the PDE in Fig. 2A) can be 
determined by differentiation of p with respect to (j). Derivative of p with respect to (j) is 
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For Z{(j)) = 1 — cos{(j)), M{(j)) would be 
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and Eq. S30 reduces to 


dp 

d(j) 

therefore 

0 = 

and 


Aw 


1 

T 


T 


COs{(j)) 


.MW 




1 


Aw 


I 

T 


T 


cos((/)*) ) 


Jo e 




MWd(/> 

Aj — 1 


Jp e~^^'Ad(j) Jo 


.MW) 


-1)- 














- 

Jo e 




MWd(j> 

— 1 




1 /Aw I 


a 


“ A-cos(/>*) /p e-MWd(j) 


Ag T T 


By using equations S34 and S28 we have 

1 e-“-^(^-i)-I 


= - 




then (j}* is 

/>*(Aw) = (</) 


(S30) 


(S31) 


(S32) 


(S33) 


(S34) 


(S35) 


= 4> PALj{4>)d4>- (S36) 

Jo 

We have used MATLAB function ”int” to plot </* versus Aw in Figure 2C using Eqs. S28, S36, and Cmax] and pmax 
versus Aw in Figure 2B by Eq. S35. 

In Eigure 3 we have plotted C{t) for two LIF neurons described by 


TmVi = Vrest - Vi + li + gijG{Vi, Vj) + VD^i^t) (S37) 

with a reset condition Vi(t'^ = !/■) if Vi(t~ > Vth)- we integrated this equation for two pulse coupled neurons and 
calculated spike count for a time window of O.bmS. Parameters are selected in agreement with biological cases as 
Tm = 20 ms, Vres = ”74 mV, Vth = —54 mV, Vr = —60 mV, I = 25.0 + A/ mV , gtj = 1 mV and D = 1.0. 


























